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We describe an empirical, self-consistent, orthogonal tight-binding model 
for zirconia, which allows for the polarizability of the anions at dipole and 
quadrupole levels and for crystal field splitting of the cation d orbitals. This 
is achieved by mixing the orbitals of different symmetry on a site with cou- 
pling coefficients driven by the Coulomb potentials up to octapole level. The 
additional forces on atoms due to the self-consistency and polarizabilities 
are exactly obtained by straightforward electrostatics, by analogy with the 
Hellmann-Feynman theorem as applied in first-principles calculations. The 
model correctly orders the zero temperature energies of all zirconia poly- 
morphs. The Zr-0 matrix elements of the Hamiltonian, which measure cova- 
lency, make a greater contribution than the polarizability to the energy dif- 
ferences between phases. Results for elastic constants of the cubic and tetrag- 
onal phases and phonon frequencies of the cubic phase are also presented and 
compared with some experimental data and first-principles calculations. We 
suggest that the model will be useful for studying finite temperature effects 

by means of molecular dynamics. 
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I. INTRODUCTION 



Solid solutions of zirconia (Zr0 2 ) containing other oxides are among the major represen- 
tatives of modern ceramic materials. The wide range of applications, including traditional 
structural refractories, fuel cells and electronic devices such as oxygen sensorsjlli testifies to 
the technological importance of zirconias. Different divalent and trivalent oxides are added 
to ZrC>2 in order to improve its thermomechanical properties, and charge-compensating va- 
cancies are thereby introduced on the anion sublattice. The macroscopic effects associated 
with the impurities are very well known,! - ! but a microscopic model which gives a theoretical 
interpretation is still missing. As a preliminary step, this paper provides a physical picture of 
the crystal thermodynamics of pure zirconia, combining the results of first principles density 
functional and semiempirical Tight Binding (TB) calculations. 

Zirconia has three zero-pressure polymorphs; these have cubic (c), tetragonal (t) and 
monoclinic (m) symmetry. The high temperature c phaselQ (Fm3m) is stable between 2570 
K and the melting temperature of 2980 K. The t structure!'! (P^/nmc), which is stable 
between 1400 and 2570 K, is closely related to the c one: the internal degree of freedom 
S shifts the oxygen ions away from the centrosymmetric positions along the X 2 mode of 
vibration (Figure [I]) and forces the c/a ratio of the unit cell to adjust. Below 1400 K the 
low-symmetry m phase0Hll (P2i/c) is thermodynamically stable. 

Besides its technological implications, the relationship between these structures is of 
fundamental interest. The mechanisms of the phase transformations, the effects of impurities 
and vacancies on them, and their relationship to the nature of the bonding still require 
explanation, and this may shed light on the properties of other, more complex oxides. 

The crystal structure of purely ionic bonded materials can be determined on the basis 
of radius-ratio rulesjl! based purely on electrostatic arguments. Because of the small size of 
the Zr 4+ ions, these rules place Zr02 on the border between the 8-fold coordinated fluorite 
structure and the 6-fold coordinated rutile one {PAz/mnni). The radius-ratio is too blunt a 
tool to account for the absolute stability of the unique 7-fold coordinated m structure. 
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The classical empirical models of zirconia are based on the a priori assumption of its 
ionicity. Empirical approaches like the Shell Model (SM) or the Rigid Ion Model (RIM) 
described the structural,0 dynamically^ and transportEHll properties of the phases on 
which they were parameterized, but failed to predict the absolute stability of the m structure. 
The most detailed of such models was developed by Wilson et aZ.,@ whose environment- 
dependent Compressible and Polarizable Ion Model (CIM-DQ) demonstrated the importance 
of the anion polarizabilities at both dipole and quadrupole levels on the energetics of zirconia. 
However, further calculations^ carried out with this model revealed that even though it 
predicted the correct energy ordering of the c, t and m phases, it predicted that the rutile 
structure should be even more stable, and this phase is never observed experimentally in 
zirconia. 

The experience gained with the CIM-DQ model suggests that a successful empirical 
model of zirconia should describe the effects of the atomic polarization, but should also go 
beyond a purely ionic description of the bonding. The partial covalent character of zirconia 
has already been postulated^ and is evident from electronic structure calculations based 
on density functional theory. In this paper we further investigate the recently proposed 
polarizable self-consistent tight binding (SC-TB) mode&~il which combines the physical 
concepts of covalency, ionicity and polarizability. Using the SC-TB model we are drawn to 
the conclusion that the covalent character of the Zr-0 bond makes a significant contribution 
to the relative energetics of different structures, which would explain the limited predictive 
power of the previous ionic models. 

There have been several previous approaches to analyzing the structural and electronic 
properties of zirconia. Boyer and Klein@ used the APW method to derive pair potentials 
with which to investigate the equation of state of the c phase. Cohen et a/.0 calculated the 
relative energetics and the elasticity using the Potential Induced Breathing (PIB) method 
based on the Gordon-Kim approach. Zandiehnadem et a/.S studied the electronic structure 
with a first principles LCAO method. The FLAPW calculations of Jansenil predicted 
for the first time the correct energetic ordering between the c and t structures at zero 



absolute temperature, identifying the double well in the potential energy that governs their 
relative stability The double well was subsequently confirmed by ab initio Hartree-Fock 
(HF) calculations,B0 but these did not predict the stability of the m structure over the t 
one. Only the very recent Density Functional Theory (DFT) calculations&B consistently 
reproduce the relative energetics of the three zirconia polymorphs at K. 

The plan of the present paper is as follow. In Section |T| we describe the model used 
in the calculations, the inclusion of the atomic polarizability in the TB framework and the 
parameterization procedure. A preliminary account of this work has been published. @ We 
have made DFT calculations of band structures of the simple structures for this purpose, 
using a new full-potential, linear muffin tin orbital method (NFP-LMTO). The predictive 
power of the new model is tested against the DFT calculations in Section [III A| where we 
study the relative energetics of zirconia. Section |I [B| focuses on the relationship between 
the c and t structures: the Landau theory of phase transformation is used to interpret the 
results of the static calculations. In Section [TV], we explore the elastic and the vibrational 
properties of the high symmetry phases. The results are summarized in the concluding 
Section. 



II. THE TIGHT BINDING MODEL 

A. Including polarizabilities in TB 

In the TB approximation the crystal wave function can be expressed as a linear combi- 
nation of atom-centered orbitals which we denote | RX): 

|^) = £cfi|RL>. (1) 

L is a composite angular momentum index L = (£, m) of the atomic orbital centered on the 
site whose position is R, n and k are the band and k— vector indices of the single particle 
wave function. For the purpose of derivation, we express the local orbitals as a product of 
a radial function and a real spherical harmonic 
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(r|RL) = / R£ (|r-R|)F L (r-R), 



(2) 



although in our empirical TB scheme the explicit functional forms of the radial wave func- 
tions are not required. To simplify the notation, we will frequently suppress the site index 
R, in which case one can take it we are referring to an atom at the origin and r is a small 
vector in its neighborhood. 

The total Hamiltonian Ji can be expressed as a sum of two terms, H = H° + Ti' . In 
traditional Self- Consistent (SC) TB, HP contains both on-site and inter-site terms. The 
on-site terms are diagonal in L, and are often taken as Hartree-Fock term values of the 
isolated atoms. The inter-site terms are the bonding integrals. The additional part of 
the Hamiltonian, Ti', is diagonal in R and L in the traditional approach (Majewski and 
VoglUfl). It controls the charge redistribution between neighboring sites which results from 
the balance between the opposite effects due to the on-site Coulomb repulsion (Hubbard U) 
and Madelung potentials. 

What is missing in the previous model is the effect of the crystal fields on the valence 
electrons, i.e. the atomic polarizability. In a preliminary account of this work@ we indicated 
how to include the polarization effects in a SC-TB formalism by adding off-diagonal terms 
^rlrl' t° the on-site blocks of the Hamiltonian. Here we describe how we make that 
extension. 

If we assume the on-site charge distribution to be localized, then its total multipole mo- 
ment Ql has a monopole contribution from the ionic core charge and a multipole (including 
monopole) contribution from the valence charge: 



As Stone^H points out, the electronic multipole moment on a site is the expectation value 
of the operator 



Ql = Q'S L0 + QI. 



(3) 



Ql = ef t Y L (T), 



(4) 
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where e is the charge of the electron. Neglecting inter-site terms like (R'L' Orx R"L") for 



R', R" 7^ R, the definition of the on-site multipole moment is therefore: 



occ. 



% = EE ° n v c 1j> ( L ' \oi\ l ") • (5) 



L'L" nk 



By invoking equations (|2|) and (^), the last factor of Eq.([5|) can be expressed as a product 
of two quantities, the Gaunt coefficients Cl>l>>l, which dictate the selection rules, and the 
integrals A^/^, which will be new parameters of the model: 



where dQ stands for the element of solid angle sin 6 d6 d<\>. The role of the Gaunt coefficients, 
which depend on the angular part of the wave function only, is to select the term with 
symmetry L arising from the coupling of the on-site orbitals V and L" . The A parameters, 
depending on the radial part of the wave function, determine the magnitude of the coupling. 
The substitution of Eq.@ in Eq.(|) defines the multipole moment of symmetry L on the 
site R. 

Having defined the on-site multipole moments, we can calculate the fields which they 
generate on all the lattice sites. The derivation uses standard results from classical electro- 
statics. The electrostatic potential is expanded in partial waves about the site: 




(6) 




(7) 



(8) 




(9) 



L 



where, using the Poisson equation, 



Vl 



47r E E^'( r, )Qr'l', 



(10) 



and 




(11) 
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The sum over L" is restricted to the values for which I" = £ + £'; Bll' are proportional to 
the well known LMTO-ASA structure constants.^ The component of electrostatic potential 
Vl couples different orbitals on a site giving the matrix elements: 

(V I W I L") =Y. v l A«"^ C VV , L . (12) 

L 

The diagonal elements of the Hamiltonian are adjusted by using a single Hubbard U 
in the standard way, which adds a term U SN-r£ to each diagonal matrix element. The 
quantities SNm are the changes in the electronic charge projected onto a site and orbital 
compared to the input, non-self-consistent charge. We use the standard Mulliken projection. 
Finally the Schrodinger equation is solved using a self-consistent iterative procedure with 
charge mixing to obtain the coefficients c^ L and hence the multipoles. 

It is useful to step back at this point and compare the above model with the Hohenberg- 
Kohn-Sham (HKS) one, whose exchange and correlation energy functional t/ xc [n] has been 
expanded to second order in the electron density n(r):il 

occ 

^hks = J- | T s + V XC + + Vj | ^ nk ) (13) 

n,k 



+U xc [n ] - Jv o xc n dr - U R [n ] + U u 

5n Sn' dr dr. 



1 ul e 2 <) 2 U- n - 



2 JJ \ I r — r'l 5n 5n' 

tiq denotes a reference electron density, which we will consider as a superposition of spherical 
ionic charges; Ts is the kinetic energy operator of the non-interacting electron gas, V IC , V H 
and Vq are the exchange and correlation, Hartree and ionic potentials calculated at the 
reference charge uq\ 5n denotes the deviation from that reference (Sn = n — uq) and n' 
refers to the electron density at r'. U u and U n are respectively the Hartree and the ion- ion 
electrostatic energies. 

Without the last term, this is simply the Harris-Foulkes functional. It generates a non- 
self-consistent TB model in which the first term is the sum of the eigenvalues while the 
second is a sum of pair potentials. If the last term is included, the total energy must be 
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minimized iteratively, and the last term now provides the self-consistency correction to the 
Kohn-Sham Hamiltonian. 

The last line of Eq. ([13]) represent the Hartree energy of the deviation from the reference 
charge, [/ H [<5n], and the second order term of the U xc Taylor expansion. We can identify 
this term in our SC-TB model as follows: 



1 rrf e 2 5 2 U XC 



2 JJ y |r — r'| SnSn' 

= W( U5N m + QnLV R L) 



SnSn'drdr'= (14) 



^ RL 

Our total energy in the SC-TB model is therefore 

occ 
ra,k 

+W(U5N* u + Q RL V nL ) (15) 

A RL 

It can be verified that, by minimizing the above expression with respect to the expansion 
coefficients in the wave functions, we recover the Schrodinger equation with the SC-TB 
Hamiltonian. 

Calculation of the forces on the ions is very straightforward once we have the self- 
consistent wave functions and multipoles. For if an ion is moved a small distance <5R, 
there is no change in total electronic energy to first order in the <5c^. Therefore we can 
calculate the force due to the change in the first term of (|i~5"D by the conventional formulae, 
using the derivatives of the non-self-consistent Hamiltonian matrix elements (see following 
section). In calculating the forces due to the last term of (|i~5"D we can hold the multipoles 
fixed and use standard electrostatics. There is no contribution to the forces from the on-site 
energy containing U. The simple form of these results for the forces in TB is a direct analogy 
with the application of the Hellmann-Feynman theorem in DFT. 
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B. Parameterization 



Each parameter of the model has been adjusted to the results of NFP-LMTO calculations, 
details of which are specified in the previous work on zirconia.@ Our TB description of 
zirconia uses a minimal basis of atomic orbitals. The oxygen atoms are modelled with 2p 
and 3s orbitals and with a fixed core charge of +4, while on the zirconium atoms there are 
4d orbitals and a core charge of +4. The purpose of the 3s orbital on the oxygen is twofold: 
to allow an extra degree of freedom for polarization, which is otherwise restricted to charge 
transfer between its 2p orbitals, and to better reproduce the structure of the conduction 
bands. 

A repulsive Born-Mayer pair potential £7 pair has been chosen in order to reproduce the 
lattice parameter and the bulk modulus of the c phase. Only the first Zr-0 coordination 
shell has been included in this interaction. 

The Hamiltonian 7i° has been adjusted to the ab initio electronic structure of the c phase 
shown in Figure (c). We chose the Goodwin-Skinner- Pettifor§ distance dependence of the 
10 hopping integrals involved. The Hubbard U have been fixed to 1 Ry. The parameters of 
the SC-TB model are collected in Table |. 

The basis set chosen reduces the number of symmetry-allowed A parameters to 4: A spp , 
A ppd , A ddd and A ddg . The first two refer to the s and p orbitals of oxygen ions, the last two 
to the d orbitals on the zirconium. 

In the highly symmetric c structure the first non spherical terms of the potential Vl on 
the cation and anion sites have g and / symmetry respectively. The latter cannot interact 
with the oxygen orbitals, the former splits the energetic levels of the zirconium d orbitals 
and A ddg determines the magnitude of the energy splitting 6e. Cubic crystal field theoryE! 
predicts the proportionality between 6e and the radial distribution of charge <r 4 > which 
is the definition of Addg given in Eq.@. Figure (a) and (b) shows the effect of the Addg 
polarization term on the band structure of the c phase: the splitting of the d bands could 
not be captured with the SC-TB without the polarizability parameters. Reasonable values 



9 



of the Addd parameter have no significant effect on any physical properties studied here, 
therefore we set it to zero. 

Less symmetric structures are necessary to parameterize the remaining A's. In the rutile 
phase, the i = 3 component of the crystal field acting on the oxygen ions splits the p levels. 
Consequently, it contributes to the width of the 2p band: this effect is controlled by A ppd 
which we adjust to match the ab initio band structure of the rutile phase. The last term 
A spp has been chosen in order to reproduce the depth of the double well in the potential 
energy of the t structure. 

III. ENERGETICS OF BULK PHASES 
A. Energy- Volume curves 

1. Zero-pressure phases 

The predictive power of the polarizable TB model has been investigated by comparing 
its results with NFP-LMTO calculations. The Energy- Volume curves calculated with the 
two methods are shown in Figure |3]. Each energy value involved the full relaxation of all the 
degrees of freedom of the structures. 

The c and the t phases were used in the parameterization procedure, therefore there is 
automatic agreement of the two methods for these crystal structures. The true prediction of 
the model is the absolute stability of the monoclinic phase. This indicates the transferability 
of the parameters between the phases. 

The rutile phase, which is not experimentally observed, has been included in the study 
because further calculations with the CIM-DQ0@ model predicted the rutile phase to be 
more stable than the monoclinic one. Figure § shows that the SC-TB model does not suffer 
from this problem, although the relative energy of the rutile phase is less than with the 
DFT. To our knowledge, the SC-TB is the first semi-empirical model which reproduces the 
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correct ordering of these polymorphs at zero temperature, including the stability of the m 
phase. 

Table [Tl] summarizes the structural properties calculated with the NFP-LMTO method 
and with the polarizable SC-TB model, comparing them with other theoretical and experi- 
mental works. The c and m lattice parameters are referred to the 12-atoms unit cell, while 
the t ones are given in terms of the 6-atoms unit cell. A comparison of the energy differences 



between the phases of zirconia calculated with different methods is given in Table [HI . 



2. High-Pressure phases 

Under pressure, the low temperature m phase transforms to an orthorhombic struc- 
ture, known as ortho I (o/), whose crystallography is still controversial. X-ray diffraction 
analysisIl'El suggests it belongs to the Pbcm space group while neutron diffraction studies!!! 
propose the Pbca space group. We carried out the calculations using the latter structure. 
The phase transition pressure strongly depends on the state of the sample and is believed 
to be between 3 and 6 GPa.&^S A second pressure-induced phase transition is observed 
around 15 GPa,0 where the oj transforms to the orthorhombic phase termed ortho II (ou). 
The latter is isostructural to cotunnite (PbCl 2 ) and belongs to the Pnam space group.0 
The pressure increases the coordination number of the zirconium atoms from 7 to 9. 

A comprehensive first-principles study of the two orthorhombic phases has apparently 
not yet been made: Stapper et a/H studied the oj structure only, while Jomard et a/S 
focused on the oji phase. 

The atomic environment of the high pressure phases is completely different to that of the 
c and t phases used in the parameterization of the TB model, therefore these orthorhombic 
structures provides a severe benchmark for the transferability of the TB parameters. 

The energy ordering of the phases predicted by the TB model is 

jjm < jjoi <l7 t <T JC < Jjou^ 
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which is the same as we obtain by combining the results of Refs. |33] and |34j. The numerical 
values of the energy differences are summarized in Table [TTT] and compare reasonably well 
with the ab initio results. The Energy- Volume curves of the orthorhombic phases are shown 
in Figure ||: all the degrees of freedom were fully relaxed and their values are collected in 
Table [TV|. 

Although the TB model predicts the correct relative energetics of the phases, it is not 
capable of describing the subtle pressure-induced phase transformation m <-> oi. Figure £| 
shows the common-tangent between the m and the on phases. As the pressure is increased, 
the model misses the correct sequence of the phases, predicting a m <-> ojj pressure- induced 
phase transformation at 5 GPa. 

B. Cubic versus Tetragonal Phases 

1. Static calculations 

The relationship between the cubic and the tetragonal phases is governed by a volume 
dependent double well in the potential energy. Since the FLAPW calculation of JanseniH 
who predicted it first, the double well has been confirmed by several other ab initio calcula- 
tions and it is now well established. 

In this section we analyze the nature of the K energy surface by combining the in- 
formation gained using two very different approaches: the NFP-LMTO method and the 
polarizable TB model. The qualitative and quantitative agreement between the results of 
the two calculations, shown in the previous section, entitles us to use the physical picture 
provided by the simpler model to interpret the ab initio results. 

Starting from the c phase, the t structure can be obtained by continuously stretching the 
unit cell along the c crystallographic direction and by displacing the oxygen columns by S 
along the tetragonal axis according to the mode of vibration (Figure [1]). We calculated 
the total energy of the crystal using the two methods, for different values of (5, c/a) at 
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several volumes. 

The energy curve exhibits a single well or a double well structure depending on the 
specific volume. At small volumes, Vi, the tetragonal distortion is energetically unfavored 
and the equilibrium structure is cubic (Figure ^[). When the cubic phase is stable, there is 
no distinct metastable tetragonal phase with which to compare its energy, so the energies 
of the two phases merge. At larger volumes, V2, a structural instability appears and the c 
structure spontaneously distorts to the t one (Figure ^). 

The curvature of the energy surfaces is related to the phase transition mechanism. It 
is clear from Figure || that |ef is positive, while is negative: this suggests that the 
phase transition is driven by the S instability and that the adjustment of the c/a ratio is a 
secondary effect. The coupling between these two order parameters will be further discussed 
when we interpret the double well using Landau Theory. 

Our LDA and TB results for the depth of the double well at the t phase equilibrium 
volume, V2, are consistent with the recent LDA values of ~ 7 mRyffl This energy barrier 
for the 6-atom unit cell corresponds to a temperature of ~ 1100 K. The same result was 
obtained by Jansenil with the FLAPW method who proposed a value of ~ 1200 K. It is 
natural that these temperatures, extrapolated from the K potential energy, underestimate 
the experimental phase transition temperature of 2570 K.i The experimentally observed 
phase transition temperature can be considered as the sum of the kinetic contributions 
of all the activated eigenmodes, while the calculated energy barrier refers to the kinetic 
contribution of the A 2 ~ eigenmode only. Even though it is reasonable to expect that at the 



phase transition the soft mode in the phonon spectra (Figure [II]) will be highly weighted 
in the total density of states, the kinetic energy kT associated with all the other modes of 
vibrations will still contribute to the measured phase transition temperature. 
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2. Physical interpretation of the double well 



What causes the c «-> t symmetry breaking? The tetragonal distortion of the oxygen 
sublattice implies the following geometrical changes: (i) Two Zr-0 bond lengths get smaller 
and two get longer but the average Zr-0 distance increases, (ii) entire columns of oxygen 
atoms shift one with respect to each other (see Fig. therefore the nearest neighbor 0-0 
distances along the column remain constant while the other 4 nearest neighbor 0-0 distances 
increase, (hi) All the Zr-Zr distances remain constant. The overall increase of both the Zr- 
O and the 0-0 bond lengths is the basis of our interpretation of the double well, founded 
mainly on electrostatic arguments. 

By adjusting the various parameters describing ionicity, covalency and polarizability of 
the TB model we can select and isolate the effects that induce the double well, but before 
doing so it is instructive to understand how a simple RIM answers to the same question. It 
has been shownEil that it is possible to reproduce the double well with a RIM in which there 
are two contributions: a repulsive short ranged pairwise interaction U pa,a and a long ranged 
electrostatic term U n . 

jjRm = j2Ae- br ^ + £ — = UPah + U*, ( 16 ) 

i<j i<j 

z is the ionic charge and is the interatomic distance between the ions % and j. 

The Zr-0 bonds increase and decrease in length in a symmetric way. As a net result, 
the centrosymmetric position of the oxygen atoms is a relative maximum of the Coulomb 
energy U n . The change in the Madelung potential caused by the tetragonal distortion is 
shown in Figure |7| (a). The overall increase of the O-Zr and 0-0 distances makes the oxygen 
sites much more sensitive to the change of the Madelung potential then the zirconium ones. 
The structural instability can therefore be interpreted as an effective way of minimizing the 
electrostatic energy of the oxygen sublattice. The repulsive Zr-0 interaction counteracts 
the structural instability driven by the electrostatics, in a way which dominates at large 
displacements because of the exponential distance dependence of this repulsion. The double 
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well shape of the energy profile is due to the different functional form of these opposing 
energy terms of Eq. (|T6D . This argument clearly depends on the strength of the repulsion, 
and does not work if the repulsion is too weak. 

It can be noticed that analogous terms are present in the TB model and a similar in- 
terpretation is tempting. However, we now have the additional effects due to polarization, 
covalency, and charge redistribution. Figure [7| (b) shows that the absolute value of the 
self-consistent equilibrium charge Q decreases on both species. Consequently, in this ap- 
proximation, the on-site energy 

lj2[ U6N li + QnLV RL \ , (17) 

A RL 

plotted in Figure [7| (c), decreases not only because of the previous geometric arguments but 
also because the charge redistribution reduces the ionic charges and therefore both the 0-0 
and Zr-Zr electrostatic interactions. 

It is interesting to note that, on the oxygen atoms, the self-consistent charge \Q e \ de- 
creases with 5 even though the total on-site potential [the sum of the Hubbard and electro- 
static terms as in Eq.flTTDI increases. This non- intuitive behavior of the charge transfer is 
due to covalency. The charge transfer is controlled both by the on-site potential and by the 
bonding integrals, which depend on the Zr-0 distance. For 5^0, the overall increase in 
the Zr-0 distance results in a decrease in the magnitude of the hopping integrals, and this 
overcomes the opposing effect of change in the on-site potential, pushing back some electrons 
from the oxygen to the zirconium sites. 

In the CIM-DQ, it was the quadrupole polarization of the O ions which stabilized the 
tetragonal structure, so it is of interest to see if it is also the development of a quadrupole 
moment in the tetragonal phase which stabilizes it within the SC-TB model. 

In fact it turns out that covalency is the main effect, although polarizability is still 
significant. The t structure is stable with respect to the c one even with a non polarizable 
SC-TB model [Figure ^ (a)]: the small energy difference is due to both ionicity and covalency 
of the crystal. The addition of the oxygen polarizability enhances the energy difference 
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between the two phases deepening and broadening the double well [Figure § (b)]. 

We can be more specific about the nature of the polarization. In the c structure, the 
first non-zero components of the electrostatic potential are Vo and V3. The latter could, 
in principle, induce an octapole moment Q3 on the anions. We truncated the multipolar 
expansion of the atomic multipole moments to the quadrupoles Q2 therefore, within this 
approximation, the ions in the c structure are not polarized. Higher order terms can be 
included in the expansion, but the overall agreement of the results with both experiments and 
first-principle calculations demonstrates that the model is already capturing the important 
physics of the system. 

As the anion sublattice is distorted, the symmetry lowering induces the £ = 1 and £ = 2 
components of the potential which couple the s and p oxygen atomic orbitals. The magnitude 
of the coupling, and therefore of the multipole moments, is controlled by the parameters 
A spp and A pp d. The latter, fixed in order to reproduce the electronic structure of the rutile 
phase, produces very weak quadrupole moments, whose contribution to the double well is 
negligible. The former controls the size of the dipole moments whose symmetric distribution 
further minimize the electrostatic energy [Figure (d)]. The total effect on the double well 
is shown in Figure ^|. 



3. Landau theory 

The c t phase transition can be interpreted in terms of the Landau Theory.il In a 
subsequent paper we plan to explore the free energy surface at T > with this formalism, 
so it is convenient to introduce it here to discuss the T = results. Experimentally, the 
mechanism of this phase transition has been very controversial and a clear description is 
still missing.&~Eil 

Theoretically, ChanEl suggested that a partial softening of an elastic constant is the 
driving force of this phase transition and, after symmetry considerations based on the elastic 
strains only, concluded that the phase transition must be of first order. We show here that 
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the inclusion of the order parameter S gives a second order phase transition. A similar 
discussion has been given by Ishibashi and Dvofak.ll 

According to the Landau Theory, the appropriate thermodynamic potential which de- 
scribes the relationship between the two phases of interest, is expanded in a Taylor series in 
one or more order parameters, in which the expansion coefficients are temperature depen- 
dent. The order parameters are non-zero in the low symmetry phase and vanish in the high 
symmetry one, providing therefore a unique way to differentiate the two phases. The terms 
involved in the Taylor expansion are invariants under the symmetry operations of the high 
symmetry phase and can be identified using group theory. 

In the case of zirconia, the c structure is unstable along the three crystallographic direc- 
tions, therefore the distortions along x, y and z have to be explicitly treated in the energy 
expansion. This suggests the following 9 order parameters, defined in terms of the strain 
tensor e, and grouped into 4 symmetry-adapted bases which spans the corresponding irre- 
ducible representations: 

S x ,S y ,S z Ti 

^xx £yy ^zz ^-1 
zz ^xx ^yy) i (^€ xx ^-yy) E 
^xyi ^yzi ^zx ^2 

A complete analysis involving all the order parameters will be done in a separate paper, 
here we simplify the total energy expansion selecting one of the three possible directions of 
the tetragonal axis. Under this hypothesis three order parameters are necessary to describe 
the c <-> t phase transition of zirconia: 5, r\ and t]q. The high temperature c phase has 
the full cubic symmetry m3m and the only degree of freedom is the hydrostatic strain 
Vo — e xx + tyy + e zz- The low-symmetry t phase is defined by the distortion of the anionic 
sublattice 5, which we define as the amplitude of the mode of vibration, and by the 
tetragonal strain r\ = (2e 22 — e xx — e yy ). 

The three order parameters can be hierarchically classified according to the amount 
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of symmetry breaking that they involve. The hydrostatic strain 770 preserves the cubic 
symmetry of the crystal. The tetragonal strain rj maintains the number of atoms in the 
primitive cell and lowers the symmetry to the point group A/mmm which still has the 
mirror symmetry operation perpendicular to the tetragonal axis. The tetragonal distortion 
S breaks this symmetry operation and involves cell doubling. Therefore, according to Landau 
theory, 5 is the primary order parameter, 77 is the secondary and 770 is the tertiary one. 

The potential energy is expanded as a power series in these order parameters around the 
equilibrium volume of the cubic phase Vq (Figure H): 



^-^ = ^ 2 + ^ 4 + M 2 » + M 2 ^+ (18) 

The elastic constants c and C\ are proportional respectively to the bulk modulus and to 
C = \{cu — C12) in the c phase described in the next section. The third order term S 3 is 
forbidden by symmetry, therefore this transition is of second order if a?, goes negative. 

The volume dependence of the order parameters can be studied by setting to zero V^U 
and V V0 U. Both the ab initio and TB results (Figure |J) confirm the analytic expressions: 



A S 2 

Cl 



Vo = ~ b fJ 2 



5 OC y/Tfr 

(19) 



V oc Vo 



These expressions show that the second-order strain terms of Eq.flrS]) are already pro- 
portional to 5 A and therefore, within the chosen order of approximation, it is not necessary 
to include third-order terms in e^-. Moreover, from the static results it is clear that the 
description of the high temperature stability of the c phase must go beyond the quasi- 
harmonic approximation. The higher the temperature, the larger the volume and, according 
to Figure |9|, the larger S and rj. Therefore, in a simple quasi-harmonic picture, a higher tem- 
perature seems to favor the t phase with respect to c, in contradiction to the experimental 
observation. 
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The parameters c\ and cq are known from the elastic properties of the crystal and have 
been calculated independently (see next section). The coefficients a 2 and have been fitted 
to the double well of an undistorted stress-free cubic crystal (in the sense rj = and rjo — 0) . 
In a similar way, b\ and 60 have been fitted to the double well of a tetragonal crystal at Vq 
(770 = 0, rj 7^ 0) and of a cubic crystal near V (77 = 0, r] 7^ 0) respectively. Figure |TD| (a) 
shows the three curves used for the fitting procedure. The agreement is very good even far 
away from the reference volume of the energy expansion [Figure |10| (b)]. This demonstrate 
that the fourth order truncation in Eq.(|l|) is sufficient to capture all the essential features 
of the K energy surface. 

Nardelli et a/.0@ have shown the crucial role played by the coupling between different 
order parameters and how it can affect the correct interpretation of the phase transformation. 
To see this we substitute the relationships (|T9|) back in Eq. (|18D : 

Cl4 6q b\ 



S' + O (5 b . (20) 



4 2 c 2ci_ 

/ b 2 b 2 \ 

The above equation shows that the coupling term ( ^ + 2cT ) can renormalize the fourth 
order coefficient, and could make it negative. In that case it would be necessary to truncate 
Eq. ( |l8l) at the sixth order term in 6, including therefore the third-order terms in the strain. 
These would then drive the phase transition making it first order. Hi! The numerical values 
of the coefficients (Table 0) allow us to estimate the amount of the coupling. We find that 
the coupling term is ~ 20% of not big enough to affect the sign of the fourth order 
coefficient and therefore the K calculations suggest that the phase transition is displacive 
of second order. 

The temperature dependence of the elastic constants might change this description and 
the final answer will be given by high temperature MD calculations which are in progress. 



IV. DISTORTIONS 



19 



A. Elastic constants 



The elasticity of c and t zirconia has been explored with the TB model. The analysis 
involved the distortion of the crystal along high symmetry directions, the calculation of 
the total energy for different values of the distortion parameter and the fit of the results 
to a polynomial. The rigidity of the crystal with respect the particular distortion applied 
has been extracted from the quadratic coefficient of the energy series expansion. For each 
strain of the t structure, we constrained the volume to the predicted equilibrium value and 
minimized the energy with respect to the internal degrees of freedom. 

Volume conserving stretches along the high symmetry directions of the c unit cell <100> 
and <111> provide C = |(cn — Cu) and C44 respectively. Extra distortions are necessary 
when the symmetry is lower: if z is the tetragonal axis, an independent set of 5 shear moduli 
were obtained by stretching along <100>, <001>, <111>, <110> and <101>. The bulk 
moduli have been obtained by fitting the Energy- Volume curves with a Birch-Murnaghan 
Equation of StateMS 

Liu et a/Jilused the slope of the acoustic branches at small wavelength of a Zr02-Y 2 03 
(15 %) system to estimate the elastic constants of the cubic phase. Kandil et aZ.0 directly 
measured the elastic constants of Yttria Stabilized Zirconia (YSZ) single crystals: the ref- 
erence values included in Table [VT| are extrapolations to 0% impurities. To our knowledge 
there is no equivalent experimental study of the elasticity of the t phase. The most recent 
values^ are measured via a powder diffraction technique on 12% Ce-doped t zirconia. 

We compare our predictions with theoretical and experimental data in Table |V|. The 
results of two other theoretical approaches, the Hartree-Fock and the PIB ones, are very 
different. As already mentioned in the Introduction, none of these calculations predicted 
the correct relative energetics of the crystal structures. Elasticity is a property of the energy 
second derivative: a good description of the energy curves is a prerequisite for reliable elastic 
constant calculations. 

The fairly good agreement of our calculations with the experiments further indicates that 
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the SC-TB model captures the main physics of the bonding. The bulk modulus, however, is 
seriously overestimated: this may not be an intrinsic limitation of the TB model, because it 
was fit precisely to the NFP-LMTO calculation, which similarly overestimates this quantity. 



B. Phonon Spectra 



In order to test the model further, as well as to give further insight into the spontaneous 
symmetry breaking of the c phase, we studied its vibrational properties. First principle 
calculationsBE! predict an imaginary frequency at the boundary of the BZ: this reinforces 
the idea that the phase transition is displacive, and driven by the softening of an optic mode. 

Our calculations were carried out with the TB model on a 96-atom unit cell. The 
eigenvalues and eigenvectors of the possible vibrational modes in that unit cell, were found 
by diagonalising the dynamical matrix which we calculated using the direct method. The 
procedure was as follows. 

Within the harmonic approximation, the potential energy $ is expanded to second order 
in powers of the atomic displacements u: 



$ = $ 



o 













fr\ 






Up 


I, K, a 


kk' 




W 




[K'J 



(21) 



We use the notation of Maradudin et a/.H k and I label respectively the atom in the 
primitive cell and the position of the primitive cell with respect to some origin. The direct 
method consists in computing the force constants <& a p via total energy and force calculations. 
In general, the atom k in the I cell is displaced by a small amount in direction a and 
the Hellmann-Feynman forces on the other atoms are recorded. These give directly the 
quadratic terms in the total energy expansion. The force constants $ aj g can be related to 
the corresponding term of the dynamical matrix D via the usual relation: 



D 



a/3 



,-27rk-x(/) 



(22) 
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M K is the mass of the atom k and k is a point in the BZ. The crystal symmetry can 
considerably reduce the number of necessary independent calculations. 0>0 

The phonon spectra plotted along the high symmetry direction < 100 > are shown in 
Figure [11]. The main feature of the spectra is the imaginary frequency of the mode 
of vibration which corresponds to the tetragonal instability shown in Figure |5|. As already 
mentioned the tetragonal instability involves cell doubling therefore the corresponding eigen- 
vector appears at the BZ border of the c phase. The soft mode at the X point v s = 5. li is 
the natural consequence of the negative curvature of the energy surface at 5 = (Figure ||). 
Setting to zero the dipolar polarizability of the anions (A spp = 0), the X^T mode is still soft, 
v s = 0.8i, but the force constant corresponding to the instability is much smaller. This is 
consistent with Figure ^| where the same effect is studied from the energetic point of view: 
the energy curve is concave at 8 = even when the oxygens are not polarizable. 

The effect of the oxygen polarizability is evident on the T\ u IR-active mode, which 
involves the rigid displacement of the two atomic sublattices. The calculated vibration 
frequency is 7.9 THz when the anions are not polarizable and 6.3 THz when the dipolar 
degree of freedom is allowed. The closer agreement of the non-polarizable result with the 
DFT frequencies of 8.1 — 8.5 THz, together with the overestimation of the bulk modulus 
suggests that the present model could slightly overestimate both the short-range repulsion 
between closed shells of electrons, responsible for the high bulk modulus, and the long range 
polarization effects which make the T lu frequency lower than the ab initio values. The results 
might be improved with a more accurate re-parameterization but the physical interpretation 
of the ab-initio results, which is the main objective of this analysis, is unlikely to change. 

Table [VII| shows the general agreement of the TB model with other calculations and 
with the experimental data. The latter are measured by Raman spectroscopy and inelastic 
neutron scattering at high temperatures on YSZ. 

Certain non-analytical terms in the dynamical matrix have been neglected, namely those 
relating to macroscopic polarization or the Berry phase. For this reason our calculations 
cannot reproduce the LO-TO splitting of 12 THz calculated by Detraux et aL@ The non- 
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analytical terms can be approximated by knowing the Born effective charge and the dielectric 
tensor, both of which could in principle be obtained from our model. This has previously 
been done in a TB framework ,0 although not for ZrC>2, and we plan to investigate the 
effect in the future. 

V. CONCLUSIONS 

We have explored the predictive power of a polarizable SC-TB model by investigating the 
crystal stability of pure zirconia. The results of this extended TB model are in overall good 
agreement with our own ab initio (NFP-LMTO) calculations and with previous experimental 
and theoretical ab initio studies. This semiempirical model has captured the basic physics 
of the relative phase stability of zirconia with a set of parameters which are transferable 
between the crystal structures. A noteworthy improvement over all previous models is 
the absolute stability of the monoclinic structure at K with respect to the usual set of 
alternatives. This demonstrates that the model is ready to deal with more complicated 
crystalline environments such as solid solutions, high temperature distortions, or interfaces. 

The TB model predicts that the covalent character of the Zr-0 bond plays a major role 
in the energetics of zirconia, more so than the polarizability of the oxygen ions. For example, 
the double well about the cubic structure, absent in a rigid ion model, exists when covalency 
is included; it is further enhanced by including also polarizability at the dipole level. We do 
not believe that the separation between covalent effects and polarizability effects is unique, 
since it depends on the choice of basis functions. Quite possibly the previous polarizable 
ion models were capturing some effects of charge redistribution which could alternatively 
be described by covalency. It remains to be seen if a model for zirconia without explicit 
covalency could satisfactorily reproduce all the structural energies. 

The Landau Theory, used to interpret the TB and ab initio results, together with the 
lattice dynamic analysis, shows that the c t phase transition is displacive of the second 
order and is driven by the softening of the X 2 ~ mode of vibration. If it had been driven by 
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a softening of the corresponding elastic constant c\\ — Cu it would have been a first order 
transition. The partial softening of the elastic constants due to the temperature could also 
in principle change the character of the phase transition. We are currently applying the 
molecular dynamics technique to understand the high temperature thermodynamic stability 
of the c phase and to explore the character of the phase transition. To this end we can 
use thermodynamic integration to go beyond the harmonic approximation. The preliminary 
results of these calculations will appear in the near future.0 

Since the valence electrons are treated explicitly within the SC-TB model we also hope 
to be able to study the effects of point defects. This would be more difficult with a classical 
polarizable ion model because of the problems associated with charge conservation and 
redistribution. 
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TABLE I. Parameters of the polarizable SC-TB model. Energy in Ry and lengths in atomic 



units. 



On site parameters 



n° s = 0.35 



U s = 1 



Hp =-0.70 

n° d =-o.io 



u d = i 



u p = 1 



Bond integrals 
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Vu> n n c d r c 
ssa -0.060 2 4.90 6.24 
spa 0.070 2 4.90 6.24 
ppa 0.050 3 4 4.90 6.24 
pp-rr -0.008 3 4 4.90 6.24 
sda -0.050 3 4.24 4.90 



V w n n c d r c 
pda -0.100 4 4.24 4.90 
pdir 0.058 4 4.24 4.90 
dda -0.050 5 6.02 6.93 
ddir 0.033 5 6.02 6.93 
ddS 0.008 5 6.02 6.93 



^spp — 0.73 



A ppd = 1.89 



Polarization terms 



A 



ddd 



Addg = 63.5 



Pair potential 
U(r) = AeH') 



A = 181.972 



1.652 
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TABLE II. Equilibrium structural parameters for the O-pressure phases of ZrC>2- The lattice 
parameters a, b, c (a.u.), and the volumes (a.u./Zr02) of the c, t, and m structure are referred 
to the 12-atoms, 6-atoms, and 12-atoms unit cells respectively. 5 denotes the internal degree of 
freedom of the t phase (see Fig. p]), (3 is the angle of the m cell in degrees, and x, y, z are the 
fractional coordinates of the non-equivalent sites in the m structure. 



Expt. a 
Refs. || 



SC-TB 
this work 



NFP-LMTO 
this work 



PW-PP 
Ref. |32 



PW-PP 



Ref. 33 



FLAPW 
Ref. |29 









Cubic 








Volume 


222.50 


213.40 


210.33 


215.29 


220.84 


217.81 


a 


9.619 


9.486 


9.442 


9.514 


9.595 


9.551 








Tetragonal 








Volume 


229.93 


217.73 


215.16 


218.69 


225.31 


218.84 


a 


6.748 


6.709 


6.695 


6.734 


6.797 


6.747 


c/a 


1.451 


1.442 


1.434 


1.432 


1.435 


1.425 


S/c 


0.057 b 


0.047 


0.051 


0.042 


0.042 


0.029 








Monoclinic 








Volume 


237.67 


222.89 


226.13 


230.51 


236.46 




a 


9.733 


9.592 


9.417 


9.611 


9.733 c 




b/a 


1.012 


1.001 


1.036 


1.024 


1.012 c 




c/a 


1.032 


1.019 


1.057 


1.028 


1.032 c 




P 


99.23 


98.00 


98.57 


99.21 


99.23° 






0.275 


0.272 


0.274 


0.278 


0.277 




VZv 


0.040 


0.027 


0.040 


0.042 


0.043 




ZZr 


0.208 


0.217 


0.212 


0.210 


0.210 




xo 1 


0.070 


0.078 


0.069 


0.077 


0.064 




yo x 


0.332 


0.336 


0.339 


0.349 


0.324 
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ZOi 

x 02 

yo 2 
zo 2 



0.345 

0.450 
0.757 
0.479 



0.342 
0.452 
0.752 
0.472 



0.338 
0.448 
0.753 
0.478 



0.331 
0.447 
0.759 
0.483 



0.352 
0.450 
0.756 
0.479 



a The experimental values of the cubic and tetragonal structures have been extrapolated at K 
using the thermal expansion data from Ref. |6| 
b At 1568 K 

c Fixed to the experimental values of Ref. || 

TABLE III. Energy differences (mRy/Zr02) between the zirconia polymorphs and the c phase 
calculated at the minimized structural parameters of Tables |l| and [V. The experimental values 
are derived from enthalpy differences at the phase transition temperature. 



AU 



t-c 



AU r 



AU°'- C 



AU°"- C 



Expt. 

SC-TB 

NFP 

PW-PP 

PW-PP 

PW-PP 



Ref. 79 



Ref. |32 
Ref. H 
Ref. |34 



-4.2 
-3.0 
-3.6 
-3.3 
-3.5 
-1.5 a -5.9 b 



-8.8 
-7.4 
-7.7 
-7.5 
-8.2 
-4.4 a -13.9 b 



-3.6 



2.8 



-5.3 



0.7 a 7.3 b 



a LDA calculation 
b Perdew-Wang GGC calculation 
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TABLE IV. External and internal degrees of freedom of the orthorhombic structures. Lattice 
parameters a, b, c in a.u., volumes in a.u./Zr02- The fractional coordinates of the non-equivalent 
sites are denoted with x, y, and z. 



Ortho I 
SC-TB 
Ref. |80| this work 



Expt. 



PW-PP 
Ref. M 



Vol. 


228.159 


218.69 


226.7 


203.54 


196.08 


212.44 


a 


19.060 


18.737 


19.060 a 


10.558 


10.541 


10.721 


b/a 


0.522 


0.520 


0.522 a 


0.596 


0.592 


0.593 


c/a 


0.505 


0.511 


0.505 a 


1.161 


1.139 


1.163 


Zir 


0.884 


0.880 


0.884 


0.246 


0.255 


0.253 


yzr 


0.033 


0.002 


0.036 


0.250 


0.250 


0.250 




0.256 


0.256 


0.253 


0.110 


0.099 


0.111 


^Oi 


0.978 


0.978 


0.978 


0.360 


0.354 


0.360 


yo x 


0.748 


0.745 


0.739 


0.250 


0.251 


0.250 




0.495 


0.509 


0.499 


0.424 


0.421 


0.425 


xo 2 


0.791 


0.784 


0.790 


0.025 


0.022 


0.023 


yo 2 


0.371 


0.371 


0.374 


0.750 


0.749 


0.750 


zo 2 


0.131 


0.130 


0.127 


0.339 


0.338 


0.340 



Expt 
Ref. |8ll 



Ortho 1 1 
SC-TB 
this work 



PW-PP 
Ref. P 



^Fixed to the experimental values of Ref. 8C 



TABLE V. Coefficients (a.u.) of the energy Taylor expansion (18). 



a 2 = 


-0.053 


bo = 


-0.062 


CO 


= 0.621 


04 = 


0.347 


h = 


-0.152 


ci 


= 0.818 
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TABLE VI. Elastic constants (GPa) of the c and t structures. 



SC-TB 
this work 



Expt. 



PIB 



Refs. 69, 7C, 71 



Ref. 27 



HF 
Ref. EC 



DFT 
Ref. |33 



310 
175 
57 



Cubic 



194 254 
167 165 
47 61 



288 
195 
180 



222 
304 
82 



268 



190 

366 

286 

180 

80 

78 



151 

327 

264 

100 

62 

59 

64 



Tetragonal 
179 
465 
326 
83 
49 
101 
156 



197 
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TABLE VII. Phonon frequencies (THz) at the T and X points of the BZ. 



SC-TB 
this work 

6.3 
15.0 



DFT 



DFT 



Expt. 



5.li 
4.5 
5.0 
12.5 
18.1 
25.0 



Ref. |72 

8.1 
17.6 
20.1 

5.8i 
4.9 
8.9 
11.0 
17.0 
21.0 



r point 



Ref. [73 

8.5 
16.5 
19.7 



X point 



5.9z 
3.5 
11.7 
11.6 
16.0 
21.0 



Refs. ||, g |83 

9.6 
18.3 
21.1 



5.1 
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FIG. 1. Cubic and tetragonal structures of ZrC-2. Light and dark circles denote oxygen and 
zirconium atoms respectively. Arrows represent the structural instability of the oxygen sublattice 
along the mode of vibration. 




FIG. 2. Band structure of cubic zirconia. In all the panels, starting from the bottom it is possible 
to identify the oxygen 2p valence bands and the unoccupied zirconium Ad bands which are partly 



(ft) TB Ad^s = 



(h) TB # 



(c) KFr-LMTO 



o,: 
o 

_ -o,: 

£5-0.4 
u 

& -0,6 
-0.8 
- 1 





X W L TXLTK 




FIG. 3. SC-TB (top) and NFP-LMTO (bottom) Energy- Volume data for the cubic (c), tetrag 
onal (t), and monoclinic (m) phases of zirconia fitted with Murnaghan equation of states. 
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_\Q I I I I I I I I I I I 

200 205 210 215 220 225 230 235 240 245 250 




200 205 210 215 220 225 230 235 240 245 250 



Volume (a.u./Zr0 2 ) 

FIG. 4. Energy- Volume curves for the monoclinic (m) and orthorhombic (oj and ojj) phases 
calculated with the TB model. 




190 195 200 205 210 215 220 225 230 235 
Volume (a.u./Zr0 2 ) 

FIG. 5. Energy- Volume curves for the c and t structures: note the convergence at small volumes 
Vi. Vq and V 2 are the equilibrium volumes of the c and t phases respectively. 
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195 200 205 210 215 220 225 230 
Volume (a.u./Zr0 2 ) 

FIG. 6. SC-TB cohesive energy vs. tetragonal distortion 5: volume and c/a dependence, (a) 
Single well at Vi=198 a.u/Zr0 2 ; (b) Double well at V 2 =218 a.u/Zr0 2 . 
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(b) 
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-*-■ c/a=0.99 * 
-— c/a=1.00 
c/a=1.01 
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r -2 



0.05 0.1 0.15 0.2 0.25 0.1 0.2 0.3 0.4 0.5 0.6 
8 (a.u.) 8 (a.u.) 



FIG. 7. 5 dependence of: (a) Madelung potential, (b) self-consistent charge Q = Q e + Q l , (c) 



Electrostatic and Hubbard energies as in Eq. fll7), (d) the same including dipoles and quadrupoles. 
The zero of energy is the top of the double well at V 2 , total energies are in Ry/formula unit, other 
quantities in a.u. /ion. 
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FIG. 8. Double well in the TB total energy at V2 : (x) no coupling between the potential and 
the oxygen atomic orbitals; (o) with dipoles and quadrupoles on the oxygen atoms. 




0.1 0.2 0.3 0.4 0.5 0.6 
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FIG. 9. Volume dependence of the order parameters calculated with the TB model: rjo is the 
hydrostatic strain of the cubic cell from the reference volume Vo, rj is the tetragonal strain of the 
cell, and 5 (a.u.) is the tetragonal distortion of the oxygen sublattice. 




FIG. 10. SC-TB total energy versus tetragonal distortion 5. (a) Fit of the data with the Landau 

energy expansion Eq,(|i~8D; (b) transferability of the coefficients at values of hydrostatic (ryo) and 

tetragonal (rj) strains different from the reference ones. 

(a) (b) 

2 



O 

B 
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+ T1 =0 T|=0 
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□ T| =0.01 T|=0 
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8 (a.u.) 8 (a.u.) 



FIG. 11. Phonon dispersion curves of cubic zirconia in the [100] direction. Closed circles are TB 
calculations, dashed lines are guides to the eye. Note the imaginary frequency of the mode of 
vibration. 
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